Aggregation kinetics of stiff polyelectrolytes in the presence of multivalent salt 



Hossein Fazli 

Institute for Advanced Studies in Basic Sciences, Zanjan 4-5195-1159, Iran 



o 
o 

<N 

<D 
CO 



Ramin GolestanianEI 

Department of Physics and Astronomy, University of Sheffield, Sheffie 

(Dated: February 1, 2008) 



S3 7RH, UK 



Using molecular dynamics simulations, the kinetics of bundle formation for stiff polyelectrolytes 
such as act in is studied in the solution of multivalent salt. The dominant kinetic mode of aggregation 
is found to be the case of one end of one rod meeting others at right angle due to electrostatic 
interactions. The kinetic pathway to bundle formation involves a hierarchical structure of small 
clusters forming initially and then feeding into larger clusters, which is reminiscent of the flocculation 
dynamics of colloids. For the first few cluster sizes, the Smoluchowski formula for the time evolution 
of the cluster size gives a reasonable account for the results of our simulation without a single fitting 
parameter. The description using Smoluchowski formula provides evidence for the aggregation time 
scale to be controlled by diffusion, with no appreciable energy barrier to overcome. 

PACS numbers: 87.15.-v,36.20.-r,61.41.+e 
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I. INTRODUCTION 

Highly charged polyelectrolytes (PEs) are known to at- 
tract each other due to electrostatic correlations brought 
about by multivalent counterions (ions of opposite 
charge) [l[ . A ubiquitous phenomenon arising from these 
correlations is the formation of collapsed bundles of stiff 
PEs which is believed to play a significant role in 

biological processes such as cell scaffolding dynamics Q 
and motility [5]. While our current theoretical under- 
standing of the process of PE bundle formation predicts 
a macroscopic phase separation, i.e. an infinitely large 
bundle @, experiments always find finite-sized bundles 
@, @1- To explain this, it has been suggested that the 
theoretically expected phase separation may be hindered 
by kinetic barriers [7], steric effects or frustration 
of the local structure with energy penalty @, 0]. The 
phenomenon of bundle formation has also been studied 
using computer simulation, which indicated a tendency 
towards a well denned finite size [TO, EU • Recently, an 
extensive study on the thermodynamic properties of a 
condensed bundle with multivalent counterions has been 
carried out, which shows that finite bundles are stable for 
an intermediate range of electrostatic couplings whereas 
at strong enough coupling the bundle could be macro- 
scopic [12j. Multivalent counterions can also induce the 
structural collapse of single semiflexible polyelectrolytes 
[l3[ and highly charged polyelectrolyte brushes Q • 

We can understand more about this phenomenon by 
studying the kinetics of the bundle formation. The an- 
gle dependent interaction between two rods has been re- 
cently studied [15[ and it has been shown that the pre- 
ferred relative orientation of the two rods has a non- 
trivial connection with whether the overall interaction 
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is attractive or repulsive. One would then like to know 
how this complicated angle-dependent interaction affects 
the fate of the filament bundle in the course of the ag- 
gregation kinetics. 

Another question could be the dominant kinetic mode 
of aggregation. A Brownian dynamics simulation using 
a model short-ranged attraction between two rods has 
shown that the rods are most likely to meet in a cross 
configuration, presumably against a barrier as studied in 
Ref. ML and then rotate and slide to adjust into a parallel 
pair [16| . Another pathway with a lower barrier has been 
suggested in which a rod will slide parallel to an exist- 
ing bundle [17[ . While these descriptions are typically in 
terms of pairs of rods or an effective interaction of a sin- 
gle filament with an already formed bundle, it is not clear 
a priori that they are feasible in such a highly correlated 
system. In particular, the orientational dependence of 
the many-body interaction between highly charged poly- 
mers at close proximity is highly complex and can lead 
to nontrivial effects, such as spontaneous formation of 
chiral structures [18[- Recently, an experiment has been 
performed to probe the kinetics of bundle formation by 
using two different fluorescent ly labeled act in filaments, 
which suggests that actin bundles dynamically exchange 
filaments with the solution [l9[. In light of the new ex- 
perimental insight into the system, it is important to set 
up a more systematic study of the kinetics of the aggre- 
gation, so that the effect of the nontrivial many-body 
correlations during the course of the aggregation can be 
better understood. 

Here, we study the kinetics of bundle formation in the 
bulk solution of stiff polyelectrolytes in the presence of 
multivalent salt using molecular dynamics (MD) simula- 
tions. We find that the PEs undergo an aggregation pro- 
cess with doublets, triplets, etc forming and subsequently 
feeding into larger clusters. We observe that the initial 
stage of the kinetics leads to formation of PE bundles 
that have a clear size selection, up to 10-11 filaments for 
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FIG. 1: (color online). Time lapse snapshots of the system 
with N p = 4 and C3.1 = 1. +3 salt ions are shown by golden 
(light) spheres, +1 counterions by red (dark), and —1 salt ions 
by gray spheres. It can be observed that when two filaments 
or two bundles meet each other the dominant crossing angle 
is 90 degrees. In the final configuration all the PEs form a 
single bundle. 



our choice of parameters. These bundles take up all the 
smaller clusters, and are relatively much more long-lived 
than the smaller ones, while larger clusters do not seem to 
appear even when there are a number of these long-lived 
filament bundles available in the solution for possible ag- 
gregation. We find that the time-dependent size distri- 
bution of the aggregates follows a Smoluchowski fioccu- 
lation kinetics. For the first three cluster sizes, namely 
single filaments, doublets, and triplets, the time evolu- 
tion of the number of such clusters shows a reasonable 
quantitative agreement with the Smoluchowski formula 
when the energy barrier is set to zero. This result shows 
that the many-body energy barrier for the formation of 
these clusters is relatively small, and that the time scale 
for the evolution is set by the diffusion of filaments in the 
course of the aggregation process. We also monitor the 
kinetics of bundle growth and find that the dominant 
mode is due to one end of a filament /bundle meeting 
another filament /bundle either in the middle mostly at 
right angle (in the form of "H") or at one end (in the 
form of ") " ) before rotating and sliding into a parallel 
packing (see Figs. [I] and [2]). We show that these modes 
can be understood from energetic considerations by cal- 
culating the angle dependence of the potential of mean 
force between two rods. 

The rest of the paper is organized as follows. Sec- 
tion [III introduces the model and details of the simula- 
tion technique used in this work. The main results from 
the simulation are presented in Sec. HIH while Sec. [TV] 
is dedicated to the potential of mean force of two rods 
at various angles and separations. Finally, Sec. [V] closes 
the paper with some discussions and remarks. 




FIG. 2: (color online). Time lapse snapshots of the system 
with N p = 16 and C3.1 = 1. The largest bundle contains 11 
PE rods (see part f of the figure) and two smaller bundles 
containing 2 and 3 PEs finally join up and make a bundle of 
5 rods (final configuration is not shown in this figure). 



II. THE SIMULATION 

In our simulations of the bulk solution, which are 
performed with the MD simulation package ESPResSo 
v.1.8 [20], N p stiff PEs are considered each composed 
of N m = 21 spherical charged monomers of charge — e 
(electronic charge) and diameter a which is introduced 
via the Lennard- Jones potential (see below). Monomers 
of each PE are bonded to each other with separation be- 
tween them being fixed at 1.1a via a finite extensible 
nonlinear elastic (FENE) potential [2l[ and the bending 
rigidity of PE chains is modeled with a bond angle po- 
tential U<f> = k(f>(l — cos</>) with kff) = 400 k^T in which 
(j) is the angle between two successive bond vectors along 
the PE chain. We use N c = N p x N m monovalent coun- 
terions with charge +e to neutralize the PEs. We also 
consider trivalent salt with N+ positive ions with charge 
+3e and N~ = 37V+ negative ions with charge — e. In 
addition to long-ranged Coulomb interaction we include 
short-ranged Lennard- Jones repulsion between particles, 
which introduces an energy scale e into the system. MD 
time step in our simulations is r = O.Olro, in which 
ro = \Jmn 2 je is the MD time scale and m is the mass of 
the particles. 

The temperature is fixed at k^T = 1.2e using a 
Langevin thermostat. We use the particle-particle 
particle-mesh (PPPM) method to apply periodic bound- 
ary conditions for long-ranged Coulomb interaction in the 
system. The strength of the electrostatic interaction en- 
ergy relative to the thermal energy can be quantified us- 
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ing the Bjerrum length £b = £/ f sT , where e is the dielec- 
tric constant of the solvent and in our simulations we have 
used it to fix the value for a via £b = 3.2a. Following Ref. 
[IH, we define the salt concentration as = N~ /N c , 
and use values in the range of cs-.i = 0.5 — 1.2. The Debye 
length Ad = 6a for c%-i = 1. 
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FIG. 3: (color online). Nk/N p versus rescaled time t/t p for (a) k = 1, (b) /c = 2, and (c) /c = 3, with W = 1. Data and error 
bars are obtained from averaging over sets of data such as those presented in Table U corresponding to various simulations with 
identical and differing values of N p . 



In the beginning of the simulations, we fix the PE rods 
in space parallel to each other with their centers arranged 
on a square lattice of spacing in the range of 10<j — 15a in 
the middle of the simulation box and all of the ions are 
free to move for 100,000 MD time steps. We then release 
the PEs and after equilibration of the system we study 
the bundle formation kinetics for 8,000,000 MD steps. 



III. RESULTS 

We do simulations for different values of the number 
of PEs (N p = 4, 9, 16, 25, 64) and find that in simulations 
with 4 and 9 PE rods, in the final configuration of the 
system a single bundle of parallel PEs forms containing 
all of them. Figure [1] shows snapshots of the system with 
N p = 4 PE rods at different times. It can be seen that 
PEs that are going to be added to a bundle of parallel 
PEs, first meet the bundle at right angle and then the 
crossing angle vanishes. 

Simulations with larger N p show the same kinetic pro- 
cess (see Fig. [2j) as well as an aggregation with small clus- 
ters forming initially and then feeding into larger clusters. 
Table [J shows the numbers of PE bundles of different 
sizes at different times for the simulation with N p = 25. 
In this table, N^ is the number of bundles containing k 
PEs and for a system with N p PEs it is subject to the 
normalization ^ k kN^ = N p . The largest bundle at long 
times contains 11 PEs in all the simulations. In the case 
of N p = 16, two bundles containing 5 and 11 rods remain 
at long times with no affinity towards each other when 
we ran the simulation for longer times. Simulations with 
larger numbers of PEs, i.e. N p = 25 (see Table [I]) and 
64 confirm that the growth of the bundles tend to be cut 
off, or slowed down beyond the time span of our simula- 
tion, when there are 10-11 filaments in the bundle. While 
at the end of the simulation there are more bundles as 
the number of rods increases, the largest bundle at long 
times remains to be made of maximum 11 rods. 



The evolution of the clusters and their distribution — as 
exemplified in Table [J for 25 rods — resembles the aggre- 
gation kinetics of colloidal particles (22[. For the case 
of spherical colloidal particles with short-ranged interac- 
tions, Smoluchowski suggested an (approximate) expres- 
sion for the number of clusters of size k in time t as |22l 
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In this equation, the characteristic time is defined as 
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where rj is the viscosity of water, no is the initial number 
density of the particles, and W — e Ea ^ kBT is an activa- 
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tion factor (up to a numerical coefficient). The values 
of Nk/Np for k = 1, 2 and 3 obtained from our simu- 
lations are plotted in Fig. [3] as a function of t/t p with 
W = 1. In these plots, we have taken r] c± 2.4y / me/<j 
(from Ref. [23]), /c B T = 1.2e, n = 0.001 a" 3 , and 
r = O.Olro = 0.0lvW 2 /e, which yields t p = 2xl0 5 rW. 
The results are compared in Fig. [3] with the Smolu- 
chowski formula of Eq. (pQ) and an agreement in the range 
of the error bars can be seen. Naturally, the fact that the 
clusters do not exceed the maximum size of 10-11 means 
that the analogy is limited to the smaller sizes. Never- 
theless, it is remarkable that the same value for t p gives a 
good agreement for N\/N p , A^/TVp and Ns/N p: without 
even a single fitting parameter being used. 



The fact that W = 1 seems to provide a reason- 
able agreement suggests that the PE rods do not ex- 
perience substantial energy barriers in their "optimal" 
kinetic paths in the course of the aggregation, although 
an exact knowledge of the numerical prefactor for W in 
the rod geometry is required for a precise determination 
of the "activation energy" E a (which may even turn out 
to be negative). Note that the low optimal energy barrier 
does not mean that energetics does not play a role in the 
kinetics, as the dominant kinetic mode is clearly a result 
of strong electrostatic interactions. 



The problem of quantifying the kinetic barrier of the 
aggregation process is a complicated one, as it is not clear 
how one can approach it better than just looking at pair- 
potentials. Here we propose that monitoring the time 
evolution of the cluster sizes could be a good alterna- 
tive for capturing the many-body essence of the aggrega- 
tion process. The typical time scale t p that is involved 
in the evolution of the cluster sizes [see Eq. ([2j)] has 
a diffusion-controlled component r]/(nok^,T) ~ £ 3 /(LD), 

where I = n 1 ^ 3 is the initial average distance between 
the filaments in the solution, L is the length of the fil- 
aments, and D is a filament diffusion coefficient. The 
second component of t p comes from the many-body en- 
ergetics of the system. A comparison of the simulation 
results for the cluster sizes with the Smoluchowski for- 
mulas in Fig. [3] suggests that the diffusion part is domi- 
nant during the early stages of the aggregation until the 
finite-sized bundles are formed. The fact that the Smolu- 
chowski plots for zero energy barrier are close to the sim- 
ulation data points, shows that the energy barriers have 
to be small. We chose to plot the zero energy barrier form 
instead of trying to fit the data to the Smoluchowski for- 
mula and deduce an energy barrier, because we do not 
know the right prefactors for the different geometry of 
rods (rather than spheres in the original Smoluchowski 
solution). However, these prefactors must not differ too 
much from unity and therefore the argument that the 
barrier must be small at this stage will be justified from 
the agreement in Fig. [3j 
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FIG. 4: (color online). Potential of mean force per monomer 
between two perpendicular rods as a function of the sepa- 
ration between their centers zi. The size of the circles on 
C3:i = 1.0 curve corresponds to the error bars. Inset: poten- 
tial of mean force between two rods which are perpendicular 
as a function of z/2 with C3.1 = 0.7. 



IV. ANGULAR DEPENDENCE OF THE 
POTENTIAL OF MEAN FORCE 

To understand the dominant kinetic mode of PE aggre- 
gation we study the different configurations of two rods 
and keep them fixed during simulation while allowing the 
counter-ions and salt ions to move freely. In these simu- 
lations, we calculate the average force on each monomer 
and obtain the total force and torque for each PE rod. By 
integrating the average force (torque) over different sep- 
arations (orientation angles), we obtain the potential of 
mean force over a displacement (rotation about an axis) 
[IH . Let us denote the coordinates of the centers of PEs 
1 and 2 as (#1,2/1,21) and (#2, 2/2? 22) respectively (see 
the schematic configuration of the rods in Fig. [4|). 

We first show that in the range of salt concentration 
we use in our simulations, there is an attractive inter- 
action between two PE rods which are perpendicular to 
each other. Assume that rod 1 is fixed on the x axis with 
its center at the origin, and rod 2 is parallel to the y axis 
with the coordinates of its center being (0, 0, z<i). Figure 
[H shows the potential of mean force, which is calculated 
as the reversible work to bring rod 2 from z^ = 7.5<r to 
closer separations as a function of ^2, for two values of 
salt concentration cs-.i = 1 and 0.5. It can be seen that 
there is an attractive interaction between two perpendic- 
ular PEs. We also calculate the potential of mean force 
between two perpendicular rods when Z2 is kept fixed at 
Z2 = 1.5<r, and 2/2 is changed (see the schematic configu- 
ration in the inset of Fig. 2|) for c%-i =0.7. We find that 
the potential of mean force decreases with increasing 7/2, 
which shows that when the two rods are perpendicular 
at close separations, the best configuration is when one 
end of one rod touches the other at right angle (although 
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FIG. 5: (color online). The potential of mean force per 
monomer as a function of for c^-i — 0.7. A local mini- 
mum can be seen around — 90°, which is also shown as the 
inset of the figure. The size of the circles in the main figure 
shows the error bar. 

the potential of mean force is only slightly lower). 

We also fix rod 1 on the x axis with its end at the 
origin and fix one end of rod 2 at (0,0,1.5a) (see Fig. 
[5]), and calculate the angular dependence of the poten- 
tial of mean force by integrating the torque with respect 
to the rotation angle (around z axis). In Fig. the 
orientational dependence of the potential of mean force 
of the two rods is shown relative to the = 90° configu- 
ration. In this figure, we can see that although the global 
minimum of the potential of mean force corresponds to 
parallel configuration of the rods, a local minimum exists 
around = 90°, which is shown in detail in the inset of 
the figure. The attractive interaction between two paral- 
lel PEs in the presence of multivalent salt appears only 
at short separations [15[ . This means that for two rods at 
larger separations, it is more likely that they attract each 
other when they have larger relative angles, and one of 
the ends is meeting the other PE in the solution. More- 
over, Fig. [5] shows that although the local minimum is 
shallow compared to the parallel-configuration minimum, 
it covers a relatively wider range of 0, and thus crossing 
at "nearly" right angle has a considerably high probabil- 
ity. For obtaining each point of Figs. [H and O we have 
used 1.5 x 10 5 MD steps for equilibration of the system 
and averaging is done over 3 x 10 5 MD steps. 



V. DISCUSSION 

In a recent experiment, the kinetics of bundle forma- 
tion of actin in the presence of multivalent salt has been 
monitored, using fluorescence microscopy, and it has been 
found that the initial stage of bundle growth tends to 
saturate when the bundles reach the size of 10-20 fila- 
ments [l9|| . The experiment suggests that the bundles 
at this stage very actively exchange single filaments with 
the solution, and do not seem to be trapped in a non- 



equilibrium state that is hindered by kinetic barriers. It 
also shows that at later stages the bundles grow longitu- 
dinally, while still keeping the same thickness pjl . A sim- 
ilar kinetic pattern is observed in our simulations, show- 
ing an initial bundle formation that appears to suddenly 
"saturate" or drastically slow down in its growth when 
the size of 10-11 is reached. It should be noted, however, 
that in our simulation we did not probe the equilibrium 
structure of the system and therefore our results do not 
provide conclusive information on the equilibrium distri- 
bution of the bundle sizes. We also note that we did 
not observe the longitudinal growth within our simula- 
tion time scale, which is consistent with the experimental 
observation that the time scale for the diffusion-limited 
aggregation (DLA) of the already formed bundles is 2- 
3 orders of magnitude longer than the bundle formation 
time scale [19|. It is worth mentioning, however, that a 
Smoluchowski type kinetics does describe even the later 
stages of the evolution in the experiment which involves 
time scale of the order of hours, but interestingly this 
rather long time scale is indeed mostly set by the dif- 
fusion component of the aggregation time scale t p [Eq. 
d2}] and the actual energy barrier to overcome is quite 
small (of the order of k&T) [19]. This suggests that 
our proposed scheme of quantifying the kinetics of the 
growth within the framework of a flocculation dynam- 
ics which takes into account both the availability of the 
constituents (controlled by diffusion) and the energetic 
barriers to overcome can help us better understand the 
problem of finite bundle formation. In other words, just 
because a part of the process is slow in absolute time scale 
it does not necessarily mean that a large energy barrier 
is involved. We would like to point out that this delicate 
issue has not been recognized in the previous literature 
of polyelectrolyte bundles. 

One may wonder whether the simple model used here 
can properly account for the physics of charged actin pro- 
tein filaments in solution. In particular, the diameter of 
F-actin is one order of magnitude larger than the value 
we have used for our model stiff polyelectrolytes. Look- 
ing at the spatial distribution of the charges on the large 
protein surface, however, one can note that they are dis- 
tributed on narrow twisted strips with a helical pitch that 
is considerably larger than the Debye length. This means 
that the effective portions of the charge distribution on 
the actin filaments that interact with each other are in 
fact not too different from thin short rods of the same 
charge density. Since the Debye screening length is much 
smaller than the pitch, one expects that the twist struc- 
ture does not matter that much at this stage. Another 
important point is that finite-sized bundle formation has 
been observed in a variety of bio-polyelectrolytes such 
as actin, microtubule etc., each of which having very 
different detailed structures Q. The generality of the 
observed phenomenon suggests that the details of the in- 
dividual systems are probably not the key determining 
factor in the formation of finite bundles. Since all of the 
bio-polymers that make finite bundles are highly charged, 
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one is naturally led to the important question of whether 
electrostatics alone can cause this effect. This is why ob- 
servation of a tendency to form finite bundles in simple 
model polyelectrolytes like ours could provide the key to 
understanding the physical mechanisms behind this phe- 
nomenon. 

Finally, we would like to remark on a similar work 
that has been recently performed by Sayar and Holm, in 
which the thermodynamic properties of condensed bun- 
dles with multivalent counterions is studied [Hj]. They 
use a hybrid MC-MD technique, and are primarily con- 
cerned with finding the ultimate equilibrium properties 
of the system, rather than the kinetic pathways of going 
towards that equilibrium which we probe. In this sense, 
the two works are complementary as they approach a 
similar problem from different perspectives and with dif- 
ferent tools. The main finding of Ref. [12[ is that depend- 
ing on the parameters the equilibrium state of such stiff 
polyelectrolytes could be both finite bundles and macro- 
scopic condensation. Their results on potential of mean 
force can be used to deduce energy barriers, which are 
in agreement with our kinetic estimates. There are also 
slight differences in the two systems which might cause 
differences in their behaviors. For example, in Ref. [l2[ 
there are no added salt ions and the system just has 
enough counterions to neutralize. This means that en- 
tropic considerations might hamper full neutralization of 



the system and thus create additional barriers (for un- 
neutralized bundles). In our case because we have both 
salt and counterions, the monovalent ions can partake the 
entropy and therefore the multivalent counterions can 
happily reside in the bundle areas, which is perhaps a 
more faithful modeling of the actual environment of such 
systems. 

In conclusion, we have studied the kinetics of the bun- 
dle formation for stiff polyelectrolytes in multivalent salt 
solutions. The distribution of cluster sizes was found to 
follow a Smoluchowski dynamics, with no appreciable en- 
ergy barrier in the optimal kinetic paths, in contrast to 
previous suggestions [7]. We also found that the domi- 
nant kinetic mode of aggregation comes from configura- 
tions with one end of a rod meeting the other rods at 
nearly right angle, and not parallel as proposed in Ref. 
[13 • These results could hopefully shed some light on 
the controversial issue of the finite size of actin bundles. 
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